n = 100;
A = zeros(n,n);
cal_eigens = zeros(n,1);
errors = zeros(n,1);

cal_time = 0;

A(1,1) = 2; A(n,n) = 2; A(2,1) = -1; A(n-1,n) = -1;
for i = 2:(n-1)
    A(i-1,i) = -1;
    A(i+1,i) = -1;
    A(i,i) = 2;
end

eigens = sort(eig(A));
epsilon = 0.001;
delta = sqrt(2/n);

for iter = 1:100
for i = 1:n-1
for j = i+1:n
    if abs(A(i,j)) > 0
        if A(i,i) == A(j,j)
            Ja = eye(n);
            Ja(i,i)= sqrt(1/2);
            Ja(j,j)= Ja(i,i);
            Ja(i,j)= -sqrt(1/2);
        	Ja(j,i)= -Ja(i,j);
            A = Ja * A * Ja';
        else  
            Ja = eye(n);
            theta = atan( 2 * A(i,j) / ( A(i,i)-A(j,j)) ) / 2;
            Ja(i,i)= cos(theta);
            Ja(j,j)= Ja(i,i);
            Ja(i,j)= sin(theta);
        	Ja(j,i)= -Ja(i,j);
            A = Ja * A * Ja';
        end
        cal_time = cal_time + 1;
    end
end
end
for i = 1:n
    cal_eigens(i) = A(i,i);
end
cal_eigens = sort(cal_eigens);
for i = 1:n
    errors(i) = abs((eigens(i) - cal_eigens(i)) / eigens(i));
end
if max(errors) <= 0.0001
    disp(iter)
    disp(max(errors))
    break
end
delta = delta / n;
end
